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We consider fluxon dynamics in a stack of inductively coupled long Josephson junctions connected 
capacitively to a common resonant cavity at one of the boundaries. We study, through theoretical 
and numerical analysis, the possibility for the cavity to induce a transition from the energetically 
favored state of spatially separated shuttling fluxons in the different junctions to a high velocity, 
high energy state of identical fluxon modes. 



I. INTRODUCTION 



THz emission from intrinsic Josephson junctions of the 
BSCCO type has received much attention recently. Sev- 
eral experiments have been reportediii^i^, in which THz 
radiation emitted from BSCCO single crystals were ob- 
served. However in most cases the detected power is 
rather small, or the frequency is rather low, or the emit- 
ted radiation is detected indirectly on an on-chip detec- 
tor. It has also been demonstrated that BSCCO can be 
considered a Josephson junction with ac Josephson effect 
even at frequencies as high as 2 THz^. Recently a very 
convincing experiment was reported® and it has attracted 
much focus and renewed experimental efforts. 

Parallel to the experimental work there has been theo- 
retical/numerical work on fluxon dynamics in in layered 
superconductors of the BSCCO typeii^. The calculations 
demonstrate that the best way to obtain THz radiation 
is by having in-phase motion of the fluxons in the differ- 
ent layers. This poses an interesting problem since the 
in-phase state of traveling fluxons is an energetically un- 
favorable state, and two fluxons of equal polarity will con- 
sequently repel each other. However, due to the disparity 
of wave speeds for in-phase and out-of-phase solutions, it 
has been shown^"^^ - that the energetically unfavorable in- 
phase state of traveling fluxons are be stable above the 
asymptotic speed of the out-of-phase mode. It has been 
assumed that the best way to obtain that is by having 
flux flow generated by a magnetic field applied parallel to 
the a-b plane, as studied theoretically in Ref. II. How- 
ever the successful experiment in Ref. 6 was done without 
a magnetic field, and it was suggested that an internal 
cavity (based on the so-called Fiskc steps) played a ma- 
jor role in the THz generation. In this paper we study 
the fluxon modes in a stack of Josephson junctions inter- 
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FIG. 1: The identical long Josephson junctions coupled to a 
cavity at x = I. 



acting with an external cavity. We derive the conditions 
under which a large amount of current is induced in the 
cavity, and under which the external cavity may induce 
bunching of the Josephson junction fluxons. 



II. THE MODEL 

Assuming that all the junctions in the stack are identi- 
cal, the equations for a stack of long Josephson junctions 
with iV + 1 superconducting layers and N insulating lay- 
ers can be written as^ 



(1) 



where the i'th element of <p, 0*, is the gauge invariant 
phase difference across insulating layer i. The N x N 
coupling matrix, S, is given by (only non-zero elements 
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are shown) 
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with S being the coupling parameter between the layers^. 
The vector J has the components 
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+ sm ( 



(3) 



where dampng parameter a = 1/ represents dissipa- 
tion, and T] is the bias-current in the z-direction. Each 
component of J is a current in the ^-direction. 

Equations (HJ-© have been written in normahzed 
units. Space x is normahzed to the Josephson penetra- 
tion depth, Aj — yJh/2e^oJ'^d' ^ and time t is normal- 
ized to the inverse plasma frequency uj^^ = yJhcjjleJ'^, 
where \Iq is the vacuum permeability, J'^ is the critical 
current of the individual Josephson junctions, d! is the 
effective thickness of the insulating layer, and cj is the 
capacitance of the individual junctions, see Refs. and 
H for details. 

The model of the Josephson stack coupled to a series 
cavity is shown in Fig. [H where L is the cavity induc- 
tance, i? is the cavity resistance, and A^cq is the total cav- 
ity capacitance. The boundary conditions for the phases 
can be written as^^ 



and 



(4) 



N 



fc=i 

where c = Ncq/cj is the normalized capacitance and 
q is the normalized charge in the cavity. Defining = 
1/ {uo^y N Lcq) to be the normalized cavity frequency and 
and Q = sjLj (NR^cq) as the quality factor, the linear 
cavity equation becomes 



fl dq 
Q'di 



N 



(6) 



1=1 



For more details on these equations see Ref. [l^ 

Two terms are present in Eq. ([5]). The first term cou- 
ples the junctions to the cavity, by equally dividing the 
cavity current between the N junctions. The second term 
represents a direct coupling between the junction through 
the capacitors cq. In a real world situation, the junctions 
would be embedded in a resonator and couple through 
electro-magnetic radiation from the edges. The second 
term then models the radiation leaving one junction and 
ending up in another junction without being reflected by 
the cavity. This is clearly not very efficient due to the 
geometry of the stack. With an efficient cavity, the sec- 
ond term is therefore expected to be much smaller than 
the first term and may safely be neglected. This may be 



justified by numerical calculations. We thus choose to 
consider 



0^(O,t) = O and 4^^{l,t) 



N 



(7) 



as the boundary conditions for Eqs. ([I])-© and ([6]). 



III. THEORETICAL ANALYSIS 

We analyze the system in Eqs. (H])-® and (g])-© in 
the case of weak inductive coupling where 
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valid to first order in IS"]. 



A. Cavity Current 

The solution to the linear cavity equation with initial 
conditions q{Q) — and q{Q) = is 



TO_ — ma 



—''-^J2'^i,{i,t')dt' 



with m± = (^1 ± iy/iQ^ ~ Ij /(2Q). The junction 

voltage atx — l, t), thus generates the cavity charge. 
We look at the case where there is one fiuxon in each 
junction and (j>l{l,t) then becomes a voltage pulse. To 
simplify the integrations, these pulses are approximated 
by delta functions, i.e.. 



(10) 



approximating voltage pulses at i = -t- 27rn/a;*, n — 
0, 1, oj* is the fiuxon shuttling frequency in junction i 
and is the phase shift of junction i. Note, that since 
the present analysis is performed in the case of small 
15*1, all the (5- functions will have approximate the same 
amplitude, A. 

With this ansatz, the cavity current becomes 
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with = i — T* — 27rn/a;* and where H{t) is the Heaviside 
step function. Limiting the analysis to the case of a high 
Q resonator, the steady state cavity current becomes 



_ n^Ac ^ eO"' cos(^''(t-r»)+y)') 

*=i yi + eQ-' -2ee-* cos (27rfi/w') 

for i — !■ oo. The phases, (p*, are determined by 
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Only fluxons shuttling with the same frequency, uj'^ = uj 
for i — 1, ...,N will be considered. In this case, Eq. 
can be reduced to 



N 



i=l 



with 



TVa/ 1 + e - 2e 5" cos (27rri/tj) 



(14) 



(15) 



and ip^ = ip for all i. Thus, the cavity current is very 
simple when the cavity has reached a steady state. Note 
that the amplitude of the cavity current for an in-phase 
mode (t* = r-' , i,j = 1, A^) is Ne. For an anti-phase 
mode (t^ = T^+i - {-ly+^TT/uj, i = 1,...,N ~ 1) the 
amplitude is w for N even and w e for N odd. 

Using Eq. the Hamiltonian of the stack of weakly 
coupled Josephson junctions is 



H 



1 



lfl[^(^D' + ^-^osf + -f,x (16) 
(0; - S{1 - 5,^nW+^ - 5(1 - 5.^iW-^) ) dx 



with being the Kronecker delta function. Using Eqs. 
(HI)-®, (HI) and ([5]) the rate of change in energy is 



N 



5(1 - <5.,jv)<^r - 5(1 - 



To determine the amplitude of the i5-functions, we re- 
quire that in the phase-locked state the energy-exchange 



of a "collision" with the boundary is the same for both a 
fluxon solution and the J-function approximation. This 
energy-exchange is given by the time-integral of the last 
term in Eq. (fT7)) . 
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5(1 - (5i,jv)0j ^ 



dt 



where ti and t2 are taken such that they cover one colli- 
sion with the boundary. 

To model a fluxon collision with the boundary, the 
following profiles are usedi^i^i 



(t)\x,t) = 4CTHan" 



_^ ( c_ sinh ((t — T*) u'y{u/c-)/c-) 
\ u cosh {{x ~ l)^{u/c-)/c-) 



(19) 

with cr* = ±1 determining the fluxon polarity, 7(1/) = 
l/Vl — being the Lorentz factor, and the lowest char- 
acteristic velocity w 1 + 25cos(7r/(A^ -I- 1)) to first 
order in |5|. We take the same fluxon polarity in all 
junctions, thus = a for all i. 

Following Refs. [Ulil, using Eqs. ([7l). (fM)) . and (flQ l) 
in Eq. ((T8|) yields 



N 



N 

E[(l-5(2-(5.,i-(5.,jv))x 
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J2cos{uj{t' -T^) + ip) (20) 



with 



C = 47r ^ ^ ' \ ^ ^, (21) 

cosh 



and where the integration was carried out from —00 to 
00 for mathematical convenience. 

Calculation of AiJ;, for the i5-function approximation 
in Eq. pU|) gives 



A ^ 



i=l 



N 



E cos (w(r* - r^ ) -f ip) 



(22) 



Requiring AiJ^ = AiJ^ determines the amplitude of 
the (5-functions to 



C being given by Eq. (PTjl and a = ±1. 



(23) 
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B. Current- Voltage Characteristics 

The asymptotic velocity, u, present in Eq. (jl9p may be 
determined similarly to what is outlined in Ref. [3 



c- . , f 7rw7(w/c_) 
— smh 

u \ 2a;c_ 



cosh 



2c_ 



(24) 



when the fluxons are shuttling with frequency to. The 
conditions for a steady state require that the energy av- 
eraged over one period is zero, thus 



f*0 + S (]ff 



(25) 



Using Eq. p?)) . AH — gives the condition 

AHb 
2IN 



(26) 
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where / is determined fromi^ 



sinh 



47rc_ 
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2wc_ 



(27) 



The current-voltage characteristics in Eq. include 
the phase ip, such that at a given bias current the system 
can adjust this phase together with the collision times, 
r', to satisfy condition ((26|) (if possible). The phase, (p, 
is related to the fluxon shuttling frequency, uj, through 
Eqs. (|13bp and (|13bp and one may thus change the fluxon 
shuttling frequency by changing the bias current. From 
Eq. (pS)) 7] is thus obtained as a function of cj. In the 
numerical simulations in section IV we shall, inversely, 
obtain to as function of rj. 



C. Bunching 

To calculate the conditions for the cavity to induce 
bunching (in-phase motion), we consider a triangular 
fluxon configuration with one fluxon in each junction, 
modeled by r* = (— l)V/(2u). The interaction energy 
between the fluxons is first calculated by considering an 
infinite line with a lattice spacing of r, thus 



7(„/c_) 



= 4ct tan ^ e 



(28) 



The interaction energy of this configuration is 

A — -\-J—00 



+ (1 - <5,^)0^+M , (29) 
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FIG. 2: from Eqs. ([13^1 and ([130, from Eq. (I33j, and 
Fi from Eq. N ^ 2, I = 8,a ^ 0.1, Q = 100, c = 

0.02, n = 0.3, S = -0.05, and r = 4 



resulting in the well-known fluxon-fluxon force 
dHi 



Fi 



+ 5 

2n 



(30) 
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The force on the fluxons from the boundary can be cal- 
culated from dH/dr using — — (— l)'(^J/2u, vahd for 
< x < Z, resulting in 



to + - 



2^ 



dHt 
dr 

N 



dt 



(31) 
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The separation between the fluxons, r, may now be 
determined from the condition 



Fr + Fh 







(32) 



Solving Eq. (|32p for r enables one to calculate the 
current voltage characteristics, Eq. (|26p . and the cavity 
current, Eq. (O, in the steady state for a triangular 
fluxon lattice. 

For the case of only two junctions, Eq. (I3ip reads 



At 



Ff, = (1 — 5)— sin 93 sin I 

2u \ u / 



(33) 



In Fig. [5] we plot (p from Eqs. (|13bp and (|13bp as well as 
Fj from Eq. ([50)) and Fb from Eq. ([55[) as a function of 
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Lo at constant r. It is seen that when the system is above 
the resonance frequency, the force from the boundary is 
negative while the Fi is positive, thus they may balance 
each other. Below the resonance frequency the two forces 
has the same sign and the only steady state solution must 
be the one where the fluxons move in anti-phase. 

The perturbation to the current-voltage characteristic 
by the cavity is contained in the AiJ^-term in Eq. ([26|) . 
given by Eq. (|^D|) . In the case of anti-phase motion for 
two coupled junctions, this term will be zero. For three 
junctions, however, it will be non-zero. Eq. (|30p and 
pip will have the same direction for w < resulting in 
anti-phase motion below the resonance frequency. This 
suggests that only for an odd number of junctions, we 
may observe an area of negative differential resistance in 
the current-voltage characteristics. 



IV. NUMERICAL SIMULATIONS 

Numerical simulations of the full non- linear Eqs. ((T|)- 
([3]) with boundary conditions ([6]) and ([7]) has been done 
using second order finite differences for the spatial deriva- 
tives and a 5th order Runge-Kutta method with adap- 
tive step size for the temporal integration^?. The spatial 
resolution was kept at 0.01 for all considered systems. 
The initial fluxon configuration had one fluxon in each 
junction, each moving in anti-phase with the one in the 
neighboring junction(s). The system was integrated un- 
til a stabilized cavity current was obtained or 20000 time 
units had passed. The stable cavity current, g, and the 
individual voltages at x = I, 4't{l,t), was analyzed using 
interpolation and FFT— to determine the most signifi- 
cant frequency in the power spectrum which is used as 
the cavity current frequency and the fluxon shuttling fre- 
quency, bJ. 

The difference in collision times, 5t — r/u^ can be 
calculated directly from the simulation using ipKljt). It 
may also be calculated analytically using Eqs. ([24]) and 
(I32p . where the latter Eq. was solved numerically for r. 
Sometimes there were multiple solutions, r^, and we have 
chosen r = min{ri}. When no solution was found in the 
interval, we used r/u = tt/oj corresponding to anti-phase 
motion. The amplitude of the stabilized cavity current 
can be determined using a simple line search in the q 
data and compared to the amplitude of Eq. p4|) using 
the value of r obtained from Eq. (|32p. The frequency 
versus applied bias current can be compared to Eq. (pS)) . 
again using r obtained from Eq. ([5^ . 

When a steady state is found, we observe that 



^cavity — fluxon 



with n 



1 or 2, with n 



giv- 



ing a very low cavity current and therefore no significant 
difference from the unperturbed system. For high bias 
currents, the system was observed to be in the n = 2 
state and a switching to n = 1 occurred when the system 
came close to the resonance in the current voltage char- 
acteristic. Below resonance frequency, the system again 
switched to the n = 2 state. To compare the analytical 
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FIG. 3: Amplitude of q (top), current- voltage (middle), and 
r/u (bottom). iV = 2, I = 8,a = 0.1, Q = 100, c = 
0.02, Cl = 0.3, and S = —0.05. lu is the fluxon shuttling 
frequency, controlled by the bias current (rj) in experiments. 



and numerical results, we only show the numerical sim- 
ulation points where a steady state was reached for the 
n = I case. 

In Figs. [3]and|lwe have used Eqs. ifM)) . (|26)) . ((30)) . and 
([3T|l / ([33 ]) for the analytical results (shown with dashed 
lines). 

Fig. [3] shows the results on the system with N = 
2, ; = 8,a = 0.1, Q = 100, c = 0.02, n = 0.3, and 
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FIG. 4: Amplitude of q (top), current- voltage (middle), and 
r/u (bottom). iV = 2, Z = 4, a = 0.1, Q = 100, c = 
0.02, Q = 0.6, and S = —0.05. uj is the fluxon shuttling 
frequency, controlled by the bias current (rj) in experiments. 



S = —0.05. Well above the resonance frequency we get 
very little current in the cavity. As the shuttling fre- 
quency approaches the cavity frequency, the cavity cur- 
rent increases and reaches a maximum at near cavity fre- 
quency but suddently drops to near zero slightly above 
the cavity frequency. Note that all numerical results are 
only shown for frequency values larger than the resonance 
frequency. This in general agreement with the findings 
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FIG. 5: Amplitude of q (top) and r/u (bottom) as a function 
of the cavity current frequency. N = 2, I = 8,a = 0.1, Q = 
100, c = 0.02, Q, = 0.3, and S = -0.45. lo is the fluxon 
shuttling frequency, controlled by the bias current (r;) in ex- 
periments. The simulation was started at high bias current 
resulting in high cavity frequency and then the bias current 
was lowered resulting in a lowering of the cavity current fre- 
quency. At line a the system switches from a high cavity 
current frequency to a low cavity current frequency state. At 
point b the system switches from a low cavity current fre- 
quency to a high cavity current frequency state. 



in Ref. 19 where only oscillators with frequencies higher 
than the resonance frequency can be syncronised. In ad- 
dition, the forces shown in Fig. [5] are seen to be directed 
in opposite directions for frequencies larger than the reso- 
nance frequency and in the same direction for frequencies 
smaller than the resonance frequency. The fluxon-fluxon 
distance will thus decrease only if the shuttling frequency 
is larger than the resonance frequency. The correspond- 
ing current-voltage characteristic thus shows a deviation 
from the case without a cavity only above the resonance 
frequency. The fluxon separation shows a similar behav- 
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ior. Exactly at the resonance frequency, where the cavity 
current is at its maximum, the system exhibit anti-phase 
motion due to the boundary force being zero and we find 
the system to be in the n — 1 state in the numerical sim- 
ulations. Slightly above resonance frequency, the separa- 
tion has a minimum and then it increases until it reaches 
at maximum at some point and then it decreases again. 
The cavity current is thus largest slightly above resonance 
frequency, since the boundary force is zero exactly at res- 
onance frequency, resulting in anti-phase fluxon motion. 

Fig. m shows the corresponding results on the system 
with AT = 2, Z = 4, a = 0.1, Q = 100, c = 0.02, f7 = 0.6, 
and 5* — —0.05. The general behavior is the same as the 
one in Fig. [3l except for the fluxon-fluxon separation. Be- 
low and at resonance frequency, we again see anti-phase 
behavior. Slightly above the resonance frequency, we see 
that the fluxon-fluxon separation has decreased to zero, 
i.e. the system has switched to a bunched state. At 
higher fluxon shuttling frequency, the fluxon are sepa- 
rated at some distance r and at some point this distance 
become so great that we again see anti-phase motion. 

The minor discrepancies observed in Figs. [3] and [3] 
between theory and numerical experiment are primarily 
caused by the two core assumptions in the perturbation 
analysis; namely the rigid collective coordinate approx- 
imation for the fluxon, and the idealized treatment of 
the fluxon reflection at the boundaries of the junction. 
Among the approximations inherent to these assump- 
tions are omission of phonons and the change in fluxon 
dynamics during reflections. We notice, however, that 
the agreement between theory and simulations is very 
good, as can be seen in the figures. 

The rather weak force induced by the cavity on the 
fluxons can only be used to obtain bunching in the weakly 
coupled case. It is, however, essential that the fluxons do 
not move in perfect anti-phase in order to induce current 
into the cavity. The top plot of Fig. [5] the amplitude of 
the cavity current is shown for a simulation with similar 
parameters as Fig. [3] but with a much higher inductive 
coupling, S = —0.45, approaching the case of intrinsic 
junctions. The simulation was started with a high bias 
current, resulting in a high fluxon-shuttling frequency 
and the bias current was gradually lowered resulting in 
lower fluxon shuttling frequencies. Eq. (|14p gives zero 
cavity current in the case of a perfect anti-phase mode. 
In the numerical simulations, however, we do not find 
zero cavity current, but rather that the system is in the 
n = 2 state, i.e. the cavity is oscillating with twice the 
fluxon shuttling frequency. As the fluxon shuttling fre- 



quency is lowered, the cavity current increases enough 
to slightly break the anti-phase motion and the cavity 
starts to oscillate at the fluxon shuttling frequency, seen 
in Fig. [5] by following the line marked with 'a' from the 
high frequency part to the low frequency part. As the 
fluxon shuttling frequency is lowered still, the amplitude 
of the cavity current increases and the fluxon separation, 
r/u, decreases. Near the resonance the cavity current 
gets smaller and the fluxon separation start to increase 
again. At some point, the fluxons start to behave erratic, 
meaning that we can not flnd a deflnitive value of r in 
the simulations and thus we do not obtain steady state 
motion. The part of the curve near co — 0.302 where the 
value of r/u oscillates heavily shows this. At some point 
the system again switches to the anti-phase motion with 
the cavity oscillating at twice the fluxon frequency, seen 
by following the 'b'-line from the low frequency part to 
the high frequency part of the figure. We have not been 
able to determine if the erratic behaviour near uj = 0.302 
is due to a too short simulation time before we give up 
finding a steady state of if this is the 'true' behaviour of 
the system. 



V. CONCLUSION 

We have analytically calculated the cavity current and 
the current- voltage relation for N weakly inductively cou- 
pled stacked Josephson junctions coupled to a resonance 
cavity. We have shown that the cavity introduces a force 
between the natively repulsive fluxons which may be used 
to obtain bunching in the weakly coupled case. The ef- 
fect is strongest for short junction, where the boundaries 
have larger influence. Our simple analysis show overall 
good agreement with numerical simulations. In the case 
of high inductive coupling our perturbation results devi- 
ate from the simulations but the overall picture is still 
consistent with the theory. 
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